\documentclass[10pt]{beamer}
\usetheme{Madrid}

\usepackage{graphicx}
\usepackage{animate}
\usepackage{hyperref}
\usepackage{amsmath,bm,amsfonts,amssymb,enumerate,epsfig,bbm,calc,color,capt-of,multimedia,hyperref}
\usepackage{ctex}
\usepackage{fancybox}
\usepackage{times}
\usepackage{listings}
\usepackage{booktabs}
\usepackage{colortbl}

\usepackage{xcolor}
\definecolor{mygreen}{rgb}{0,0.6,0}
\definecolor{mygray}{rgb}{0.5,0.5,0.5}
\definecolor{mymauve}{rgb}{0.58,0,0.82}

\usepackage{fontspec} % 定制字体

\usepackage{verbatim}
\lstset{ %
  backgroundcolor=\color{white},      % choose the background color
  basicstyle=\footnotesize\ttfamily,  % size of fonts used for the code
  columns=fullflexible,
  tabsize=4,
  breaklines=true,               % automatic line breaking only at whitespace
  captionpos=b,                  % sets the caption-position to bottom
  commentstyle=\color{mygreen},  % comment style
  escapeinside={\%*}{*)},        % if you want to add LaTeX within your code
  keywordstyle=\color{blue},     % keyword style
  stringstyle=\color{mymauve}\ttfamily,  % string literal style
  frame=none,
  rulesepcolor=\color{red!20!green!20!blue!20},
  % identifierstyle=\color{red},
  language=c++,
}


\title{有限元方法}
% \subtitle{求解poisson方程的示例}
\author{倪金义、庄启良}
\centering
\date{\today}


\begin{document}
\thispagestyle{empty}	%使封面没有导航条
\maketitle

% ------------------目录页
\begin{frame}           %生成目录页，目录太长时加选项[shrink]
  %	\setcounter{page}{0}	%setcounter似乎对beamer无效
  \addtocounter{framenumber}{-2}%---------位置放在beginframe之后，不然无效
  \frametitle{目录}
  \thispagestyle{empty}
  \tableofcontents        % 也可以插入选项 [pausesections]
  % ----------------------列目录时，隐藏所有的小节
  % \tableofcontents[hideallsubsections]
\end{frame}





\section{有限元的基本步骤}
\begin{frame}{有限元方法的基本步骤}
  \begin{block}{}
    \begin{itemize}
    \item 生成网格：$make\_grid()$
    \item 选择单元类型、设置自由度：$setup\_system()$
    \item 构建刚度矩阵和右端项、处理边界条件：$assemble\_system()$
    \item 求解：$solve()$
    \item 结果可视化：$outputData()$
    \end{itemize}
  \end{block}


\end{frame}


\subsection{$make\_grid()$}
\begin{frame}[fragile]{生成网格：$make\_grid()$}
  \begin{block}{}
    \begin{lstlisting}
      //在test中执行
      vector<Point<2>> p = {{-1.0, -1.0}, {1.0, -1.0}, {1.0, 1.0}, {-1.0, 1.0}};
      Rectangle2D domain(p);
      vector<int> seg({std::atoi(argv[1]), std::atoi(argv[2])});
      Mesh<2>* m = new Q1Mesh(domain, seg);
    \end{lstlisting}
  \end{block}

\end{frame}

\begin{frame}[fragile]{Mesh.h}
  \begin{block}{}
    \begin{lstlisting}
      TEMPLATE
      class Mesh{
        private:
        vector<Grid<DIM>> gridList;/*< 存放所有的Grid*/
        vector<pair<int, int>> IndexofGrid ;/*< 存放grid上的所有index */
        vector<int> seg;/*< 存放各维度分段数*/
      };
    \end{lstlisting}
  \end{block}
  \begin{block}{}
    \begin{lstlisting}
      TEMPLATE
      class StructedMesh : public Mesh<DIM>{
        public:
        /**
        * @brief 给定模板单元的阶数，在每个Grid上划定自由度
        * @details 实现二维一阶、二阶模板单元
        * 
        * @param _order_of_element 模板单元的阶数
        */
        void distributeDofs(int _order_of_element);
        private:
        vector<int> seg;/*< 存放各维度分段数*/
      };
    \end{lstlisting}
  \end{block}

\end{frame}


\begin{comment}
\begin{frame}[fragile]{Mesh.h}
  \begin{block}{\href{http://web.mit.edu/easymesh\_v1.4/www/easymesh.html}{http://web.mit.edu/easymesh\_v1.4/www/easymesh.html}}
    \begin{lstlisting}
      class EasyMesh :public  UnstructedMesh<2>
      {
        private:
        string _MeshName;
        long int _ndofs;
        long int _nele;
        long int _nbndedge;
        vector<Dofs<2>> _DofsList;
        vector<Dofs<2>> _BndDofs;
        vector<vector<int>> _ElementData;
        vector<vector<int>> _BndEdges; 

        public:
        EasyMesh();
        EasyMesh(string meshname);
        void InputDofData();
        void InputEdgeData();
        void InputEleData();
      };      
    \end{lstlisting}
  \end{block}
\end{frame}
\end{comment}


\subsection{$setup\_system()$}
\begin{frame}[fragile]{选择单元类型、设置自由度：$setup\_system()$}
  \begin{block}{}
    \begin{lstlisting}
      /*Enumerating is done by using distributeDofs(),2-order element*/
      _mesh.distributeDofs(2);/*<  Q1网格对应Q1单元，Q2网格对应Q2单元 */
      _element.setMesh(_mesh);
      /**
      * set the sizes of system_matrix, right hand side vector ans the solution vector;
      */
      system_matrix = Eigen::SparseMatrix<double>(int(_mesh.getTotalNumDofs(2)),int(_mesh.getTotalNumDofs(2)));
      system_matrix.setZero();
      system_rhs = Eigen::VectorXd::Zero(int(_mesh.getTotalNumDofs(2)));
      solution = Eigen::VectorXd::Zero(int(_mesh.getTotalNumDofs(2)));
    \end{lstlisting}
  \end{block}

  
\end{frame}

\begin{frame}{构建刚度矩阵和右端项、处理边界条件：$assemble\_system()$}
  \begin{block}{poisson方程}
    \begin{align*}
      -\Delta u &= f \qquad\qquad & \text{in}\ \Omega,
      \\
      \alpha u+\beta\frac{\partial u}{\partial n} &= g \qquad\qquad & \text{on}\ \partial\Omega.
    \end{align*}

    We will solve this equation on the square, $\Omega=[−1,1]^2$, only consider the particular case $f(x) = 1,\alpha = 1,\beta = 0,g = 0$.
  \end{block}

  \begin{block}{weak form}
    The weak form is obtained by multiplying the equation by a test function $\varphi$ from the left  and integrating over the domain $\Omega$.
    \begin{align*}
      -\int_\Omega \varphi \Delta u = \int_\Omega \varphi f.
    \end{align*}

  \end{block}    
\end{frame}

\subsection{$assemble\_system( )$}
\begin{frame}{构建刚度矩阵和右端项、处理边界条件：$assemble\_system（ ）$}
  \begin{block}{分部积分、散度定理}
    \begin{align*}
      \int_\Omega \nabla\varphi \cdot \nabla u
      -
      \int_{\partial\Omega} \varphi \mathbf{n}\cdot \nabla u
      = \int_\Omega \varphi f.
    \end{align*}
    The test function $\varphi$ has to satisfy the same kind of boundary conditions, so on the boundary $\varphi = 0$.
    $$(\nabla\varphi, \nabla u) = (\varphi, f).$$
  \end{block}
  \begin{block}{离散化}
    We seek an approximation $u_h(\mathbf{x})=\sum_j U_j \varphi_j(\mathbf{x})$, where the $U_j$ are unknown expansion coefficients and $\varphi (x)$ are the finite element shape functions we will use.
    \begin{align*}
      (\nabla\varphi_i, \nabla u_h)
      = (\varphi_i, f),
      \qquad\qquad
      i=0\ldots N-1.
    \end{align*}
  \end{block}
\end{frame}

\begin{frame}{构建刚度矩阵和右端项、处理边界条件：$assemble\_system()$}
  \begin{block}{矩阵形式}
    Find a vector U so that $A U = F$,
    where the matrix A and the right hand side F are defined as
    \begin{align*}
      A_{ij} &= (\nabla\varphi_i, \nabla \varphi_j),
      \\
      F_i &= (\varphi_i, f).
    \end{align*}
  \end{block}
  \begin{block}{拼装}
  \end{block}
\end{frame}


\subsection{solve()}
\begin{frame}[fragile]{求解：$solve()$}
  \begin{block}{求解}
    \begin{lstlisting}
      Eigen::ConjugateGradient<Eigen::SparseMatrix<double>> Solver_sparse;
      Solver_sparse.setTolerance(1e-12);
      Solver_sparse.compute(system_matrix);
      solution = Solver_sparse.solve(system_rhs);
    \end{lstlisting}
  \end{block}

\end{frame}

\subsection{$outputData()$}
\begin{frame}{结果示例}
  \begin{figure}[H]
    \begin{minipage}{0.48\linewidth}
      \centerline{\includegraphics[height=4.8cm]{output_eq1.png}}
      \caption{eq1, Q2, dofs = 17689}
    \end{minipage}
    \hfill
    \begin{minipage}{0.48\linewidth}
    	\centerline{\includegraphics[height=4.8cm]{output_eq2.png}}
    	\caption{eq2, Q2, dofs = 17689}
    \end{minipage}
  \end{figure}
\end{frame}

\subsection{}
\begin{frame}{收敛阶对比}
	\begin{figure}[H]
		\begin{minipage}{0.48\linewidth}
			\centerline{\includegraphics[height=7cm]{err_q1.png}}
			\caption{eq1, Q1收敛阶}
		\end{minipage}
	\end{figure}
\end{frame}

\subsection{}
\begin{frame}{收敛阶对比}
	\begin{figure}[H]
		\begin{minipage}{0.48\linewidth}
			\centerline{\includegraphics[height=7cm]{err_q2.png}}
			\caption{eq1, Q2收敛阶}
		\end{minipage}
	\end{figure}
\end{frame}

\subsection{}
\begin{frame}{收敛阶对比}
	\begin{figure}[H]
		\begin{minipage}{0.48\linewidth}
			\centerline{\includegraphics[height=7cm]{err_eq2_q2.png}}
			\caption{eq2, Q2收敛阶}
		\end{minipage}
	\end{figure}
\end{frame}



\subsection{solve()}
\begin{frame}[fragile]{时间效率比较}	
	\begin{block}{测试结果}
		\begin{lstlisting}[basicstyle=\tiny]
Q1 修改前
Number of degrees of freedom: 66049
setup_system takes  65.41 ms
assemble_system takes  422042 ms
dealwithBoundaryCondition takes  1.496 ms
solve takes  18816.5 ms
Error l2 norm is :0.000541433, error max norm is :4.04575e-06

Q1 修改稀疏矩阵存储
Number of degrees of freedom: 66049
setup_system takes  85.507 ms
assemble_system system takes  59989.7 ms
dealwithBoundaryCondition  system takes  2.133 ms
solve system takes  22433.5 ms
Error l2 norm is :0.000541433 ,error max norm is :4.04575e-06

Q1 修改稀疏矩阵存储 修改mesh、element、femspace等数据结构
Number of degrees of freedom: 66049
setup_system takes  4.504 ms
assemble_system takes  46631.8 ms
dealwithBoundaryCondition takes  767.8 ms
solve takes  21755.2 ms
Error l2 norm is :0.000541433, error max norm is :4.04575e-06\end{lstlisting}
	\end{block}
\end{frame}




\begin{frame}[fragile]{时间效率比较}
	\begin{block}{测试结果}
		\begin{lstlisting}[basicstyle=\tiny]
Number of degrees of freedom: 263169
setup_system takes  17.882 ms
assemble_system takes  210990 ms
dealwithBoundaryCondition takes  3278.19 ms
solve takes  157124 ms
Error l2 norm is :0.000270714, error max norm is :1.01148e-06

Number of degrees of freedom: 1050625
setup_system takes  71.723 ms
assemble_system takes  758559 ms
dealwithBoundaryCondition takes  14157.3 ms
solve takes  1.236e+06 ms
Error l2 norm is :0.000135362, error max norm is :2.5297e-07\end{lstlisting}
	\end{block}
	\begin{block}{分析}
		\begin{itemize}
			\item 调整稀疏矩阵后组装时间大幅降低，调整数据结构后建立时间大幅降低（非主要）
			\item 调整数据结构后，由于边界自由度不再由mesh直接得出，改为由mesh->element->femspace，用set存储，时间增加，但依然线性且不是主要矛盾
			\item 目前时间非线性部分主存在于由矩阵求解
		\end{itemize}
	\end{block}
\end{frame}


\begin{frame}
	\huge{\centerline{The End}}
\end{frame}


\end{document}

